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Abstract 

We report on a measurement of the cosmic ray energy spectrum with the IceTop air shower array, 
the surface component of the IceCube Neutrino Observatory at the South Pole. The data used 
in this analysis were taken between June and October, 2007, with 26 surface stations operational 
at that time, corresponding to about one third of the final array. The fiducial area used in this 
analysis was 0.122 km^. The analysis investigated the energy spectrum from 1 to 100 PeV measured 
for three different zenith angle ranges between 0° and 46°. Because of the isotropy of cosmic rays 
in this energy range the spectra from all zenith angle intervals have to agree. The cosmic-ray 
energy spectrum was determined under different assumptions on the primary mass composition. 
Good agreement of spectra in the three zenith angle ranges was found for the assumption of pure 
proton and a simple two-component model. For zenith angles 6 < 30°, where the mass dependence 
is smallest, the knee in the cosmic ray energy spectrum was observed between 3.5 and 4.32 PeV, 
depending on composition assumption. Spectral indices above the knee range from —3.08 to —3.11 
depending on primary mass composition assumption. Moreover, an indication of a flattening of 
the spectrum above 22 PeV were observed. 
Keywords: cosmic rays, energy spectrum, IceCube, IceTop 



Preprint submitted to Astroparticle Physics 



February 15, 2012 



1. Introduction 

Almost 100 years after the discovery of cosmic rays, their sources and acceleration mechanisms 
still remain mostly unknown. The energy spectrum of cosmic rays as measured by various exper- 
iments follows a relatively smooth power law with spectral index 7 w —2.7 up to about 4PeV, 
where it steepens to 7 w —3.1 [T. While this feature in the spectrum called "knee" is well es- 
tablished, its origin remains controversial [2J. Most models to explain the knee involve a change 
in chemical composition of cosmic rays in the energy region above the knee. Such a change has 
been observed by various experiments [3_ but systematic uncertainties are too large to discriminate 
individual descriptions. Features in the all-particle cosmic ray energy spectrum and their chemical 
composition bear important information on the acceleration and propagation of cosmic rays. The 
measurement of the cosmic ray energy spectrum and composition is the main goal of the IceTop 
air shower array. 

IceTop is the surface component of the IceCube Neutrino Observatory at the geographic South 
Pole. Installation of IceCube and IceTop was completed at the end of 2010, with 86 IceCube strings 
and 81 IceTop stations deployed covering an area of about Ikni^ and a volume of about Ikm'^. 
IceTop was designed to measure the energy spectrum and the primary mass composition of cosmic 
ray air showers in the energy range between 5 • 10^'' eV and 10^^ eV. 

The average atmospheric depth at the South Pole is about 680g/cm^. IceTop is therefore 
located close to the shower maximum for showers in the PeV range (for vertical protons about 
550g/cm^ at 1 PeV to 720g/cm^ at lEeV). This has the advantage that local shower density 
fluctuations are smaller than at later stages of shower development. 

In this paper, we present the first analysis of IceTop data on high-energy cosmic rays and a 
measurement of the cosmic ray energy spectrum. This analysis is based on air shower data taken 
with the IceTop surface stations. The data were taken between June and October 2007 with 26 
IceTop stations operating, which comprise about 1/3 of the complete detector. 

Section |2] of this paper gives an overview over the IceTop array, and the processing and cal- 
ibration of tank signals, which are the basis for reconstructing air showers. Section [3] describes 
the dataset and run selection criteria. Section [4] introduces event reconstruction, and in Section [5] 
simulation of air showers and of the IceTop tank response are presented. In Section [6] the final 
event selection and detector performance are discussed. Section [7] describes the determination of 
the primary energy, whereas systematic uncertainties are discussed in Section [8] In Section [9] the 
results are presented and discussed. 
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Figure 1: Layout of the IceTop air shower array. Colors indicate the year of deployment and the 26 stations installed 
in 2007 are highlighted. 



2. The detector 

2. 1 . Ice Top 

The IceTop air shower array is the surface component of IceCube, covering an area of about 
1 km^ with 81 detector stations above the 86 IceCube strings. The stations are mostly located next 
to IceCube strings with a average spacing of 125 m, except for three stations placed as an infill 
with a smaller spacing in the central part of the detector, in order to lower the energy threshold 
of the detector to about 100 TeV. By 2007, 22 IceCube strings and 26 IceTop stations had been 
deployed. These stations are highlighted in Fig. [T| which shows the layout of the IceTop air shower 
array in its final configuration. 

Each station consists of two ice-filled tanks separated from each other by 10 m. The two tanks of 
each station are embedded in snow with their tops aligned with the surface in order to minimize the 



accumulation of drifting snow (see Section 2.5) and to protect the ice from temperature variations. 

The tanks are cylindrical with an inner diameter of 1.82 m, and are filled with transparent ice to 
a depth of 90 cm (see Fig. [2| . The inner tank walls are covered with a diffusely reflective zirconium 
coating. The first four stations deployed in 2005 and four tanks of the infill have a Tyvek liner with 
a higher reflectivity. This difference affects amplitude and pulse width of detected tank signals, 
since the higher reflectivity reduces Cherenkov photon absorption, leading to longer pulses. 

Each tank is equipped with two 'Digital Optical Modules' (DOMs) [1] to record Cherenkov light 
generated by charged particles passing through the tank. The DOMs are identical to those used 
in other IceCube components and consist of a 10" photomultiplier tube (PMT) [5_, plus electronic 
circuitry for signal digitization, readout, triggering, calibration, data transfer and various control 



4 




Wooden lid 



Perlite *- 58 cm 



40 cm 




•Insulation 



Diffusely reflective liner 
(Tyvek / Zirconium) 



Figure 2: Cross section of a tank showing the tank geometry with insulation and position of the DOMs. The center 
of the ice surface between the two DOMs is used as tank position by reconstruction algorithms. 

functions. The two DOMs in each tank were operated at different PMT gains, 5 • f 0^ (high-gain 
DOM) and 5 • fO^ (low-gain DOM), to enhance the dynamic range. This resulted in a linear 
dynamic range from f to more than fO'"^ photoelectrons (PE). During the data taking period used 
in this analysis all f 04 DOMs in the 26 IceTop stations were fully operational. 

2.2. Trigger and data acquisition 

A DOM records PMT signals autonomously. A signal is recorded if it surpasses a certain 
discriminator threshold, which in the case of IceTop was set to 22 mV for the high-gain DOMs 
(corresponding to about 20 pe) and 12 mV for the low-gain DOMs (corresponding to about 180 pe). 
The exact charge threshold depends on the pulse shape, which is determined by the arrival times of 
photoelectrons. After triggering, the delayed PMT pulse is sampled by 'Analog Transient Waveform 
Digitizers' (ATWDs) with three different gain channels (nominal gains are 0.25, 2, and 16) in 128 
bins with a width of 3.3 ns, corresponding to a total sampling time of about 422 ns. The analog 
samples are then digitized to 10 bits accuracy. 

Up to this point, signal recording happens independently in each DOM. To reduce the high 
trigger rates in high-gain DOMs (^2 kHz), which are mostly from low-energy showers, a hardware 
'local coincidence' between the high-gain DOMs in the two tanks of a station is required to initiate 
the readout and transmission of DOM data to the counting house (IceCube Lab). The digitizing 
process is aborted if the high-gain DOM in the neighboring tank does not also measure a signal 
above threshold within a time window of ±1 /is. The IceTop trigger condition is satisfied, if six or 
more DOMs report a (local coincident) signal within a time window of 5 /iS, which initiates readout 
of all DOMs from 10 /iS before the first until 10 /xs after the last of the six DOM triggers which 
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Muon Spectrum of D0M(61, 61) 




Figure 3: Left: A typical IceTop waveform. The blue horizontal line marks the baseline and the near vertical green 
line indicates the extrapolation of the leading edge yielding the signal time marked by the red circle. The baseline 
is below (dashed line) due to droop. Right: A typical charge spectrum receded for the VEM calibration. The 
spectra are fitted with an empirical formula to determine the peak position (see text). 

initiated the readout. The requirement of 6 DOMs means that at least two stations had to trigger. 
In 2007, the total IceTop trigger rate was about 14 Hz. 

2.3. Charge extraction and calibration 

Figure [3] shows a typical waveform measured in IceTop. While waveforms are recorded in 
three ATWD channels, this analysis used only the highest gain unsaturated (less than 1022 ADC 
counts) channel. In this analysis only the integrated charge and the signal time were used. Before 
a waveform was integrated, its baseline was subtracted by determining the average value in bins 
83 to 123 highlighted in the figure. The undershoot is caused by droop introduced by the ferrite- 
core transformer used to couple the photomultiplier tube to the DOM's front-end electronics. The 
signal time ('leading edge time') was defined by extrapolating the steepest rise of the waveform 
before the maximum down to the baseline. The absolute time scale of a DOM is calibrated with 
respect to all other DOMs to an accuracy of about 2 ns RMS |6]. 

The charge produced by a single photoelectron, the amplifier gains and the digitizers are cali- 
brated in a procedure common to all IceCube DOMs However, the signal response to a particle 
of a given type and energy traversing the tank, expressed in photoelectrons, differs from tank to 
tank, due to differences in ice quality and reflectivity of the tank walls. Therefore, the signal of 
each tank is converted to a common unit called 'Vertical Equivalent Muon' (VEM). Calibration 
was done by recording charge spectra of DOMs in dedicated calibration runs with all DOMs op- 
erated at a gain of 5 • 10^ and without requiring local coincidence (for an example see Fig. |3] 
right). These charge spectra show a clear peak due to penetrating muons above a background of 
electrons and photons. The spectra are fitted by the sum of a function describing the muon peak 
and an exponentially falling background term. Measurements with a portable scintillator telescope 
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mounted on top of tanks, restricting muons to nearly vertical angles of incidence, indicated that 
the muon peak lies about 5% lower than for the full angular range. Simulation studies confirmed 
that restricting the angles of incidence of muons shifts the peak position by about 5% [7,. The 
scaled peak is referred to as 'VEM peak'. For a given DOM the VEM unit can be expressed in 
terms of number of photoelectrons. These values average 120 and 200 photoelectrons for the low 
and high reflectivity tanks (see above), respectively. 

For the 5-month run, 15 calibration runs were used. Between two consecutive calibration runs, 



the charge calibration was assumed to be stable (see also the discussion in Section 8.4 1. 



2.4- Atmospheric conditions 

Variations of the atmosphere influence the development of air showers and thus the signals 
measured in IceTop. Since IceTop is below the shower maximum for all energies of interest in 
this analysis and for all primary masses, an increase of the atmospheric overburden leads to an 
attenuation of shower sizes. Atmospheric overburden is related to ground pressure p as Xq — p/g, 
where g — 9.87 m/s^ is the gravitational acceleration at the South Pole. While there is some annual 
variation of the ground pressure, it mostly varies on shorter time scales on the order of days. 

Besides ground pressure, the altitude proflle of the atmosphere, dXy{h)/dh, also influences the 
development of air showers. This altitude proflle has a pronounced annual cycle because the cold 
atmosphere during the winter months is much denser than the warmer atmosphere of the summer 
months. The data used in this analysis were mostly taken during the winter months. 

In the simulations used to interpret the air shower data a model of the South Pole atmosphere is 
used, which should represent the average atmosphere during the data taking period. Nevertheless, 
variations of the atmosphere around the average lead to an additional uncertainty on the measured 
energy spectrum. These systematic uncertainties will be discussed in Section |8.2[ 



2.5. Snow 

During installation, IceTop tanks are embedded in snow up to the upper surface of the tanks. 
Depending on location, surrounding surface and structures, each tank is covered by accumulated 
layers of snow of varying thickness. Each year the amount of snow on the IceTop tanks grows by 
an average of 20 cm. 

As shown in Fig. |4] the snow height for the analyzed data varied mostly between and 30 cm, 
except for four stations close to a building, which are covered by 60 to 90 cm of snow. The average 
snow height was 20.5 cm in January, 2007. 

The snow has an average density of 0.38 g/cm"^, depending on snow height and location. The 
snow on top of and around the tanks influences the response to air shower particles penetrating the 
tanks and needs to be taken into account in simulations and for the determination of the shower 
energy. 
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Figure 4: Snow heights on top of IceTop tanks measured in January 2007. All 20 newly deployed tanks had no snow 
on top, the average snow height was 20.5 cm. The dashed histogram is the snow height distribution on top of the 
same tanks measured one year later. 

3. Data set and data selection 

Event filtering and data transmission. The data used in this analysis were taken between June 1st 
and October 31st, 2007. The analysis was performed using a data sample which was transferred 
via satellite to the IceCube data center at UW Madison with limited bandwith. Due to these 
bandwidth constraints, events with less than 16 participating DOMs were prescaled by a factor 
of 5. Events with 16 or more DOMs were transmitted at a rate of 0.9 Hz and the small events at 
a rate of 2.5 Hz. 

Run selection. In order to ensure detector stability and data quality, the following criteria were 
applied to runs which were used in this analysis: 

• The run was longer than 30min. A normal detector run lasted 8 hours, and nearly all runs 
that were aborted after a short time encountered some sort of problem. 

• All DOMs were running stably. 

• After correction for atmospheric pressure variations the trigger and filter rates were stable 
and within ±5% agreement with the previous good run. Pressure correction was done by 
fitting the relation between ground pressure p and rate R with an exponential function, 
R{p) ^ exp(— /3p), yielding a barometric coefficient /3 = 0.0077/mbar Then, the rates 
were corrected to the average South Pole ground pressure of 680 mbar: 

^corrected = R exp(/3 {p - 680 mbar)) . (1) 
These cuts reduced the livetime by about 10%. 
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Event cleaning. Before starting the reconstruction, events were cleaned based on a few simple 
timing criteria. In case both DOMs of a tank triggered, the tank signal was rejected if the time 
difference between the two signals was greater than 40 ns. The analysis used only one signal per 
tank. For each high-gain DOM a saturation threshold was determined from a comparison of signals 
that triggered both DOMs in a tank. Signals with less charge were taken from the high-gain DOM. 
If the charge exceeded the saturation threshold, the charge measured by the low-gain DOM was 
used and the time was determined from the high-gain signal. Furthermore, a tank signal was also 
rejected if only the low-gain DOM triggered and the high-gain DOM was missing. 
Then, a maximum time difference of 

\tA-tB\< ^^^'^""^ + 200 ns (2) 
c 

between signals in tanks A and B of the same station was required. Here, and tg are the 
signal times in the two tanks and xa and xb are the tank locations. The tolerance of 200 ns was 
introduced in order to account for shower fluctuations. Finally, stations were grouped in clusters, 
such that any pair of stations i and j in the cluster fulflUed the condition 

|i,_t^.| < ^i_^+200ns. (3) 

The station position Xi is the center of the line connecting its two tanks, and ti is the average time 
of the tank signals. In each event, only the largest cluster of stations was kept. 

Only about 10% of events were affected by this event cleaning, and about 2.5% of events 
dropped below the threshold of 5 stations required for reconstruction. 

Charge-based retriggering. In order to reduce uncertainties due to the description of the detector 
threshold in the simulation, all events were retriggered to a common threshold based on total 
registered charge. All pulses with a charge below S'tiu- — 0.3 VEM were removed, and afterwards 



the local coincidence conditions (see Section 2.2) were re-evaluated discarding all pulses that no 



longer fulfilled this condition. This procedure was applied to both experimental and simulated 
data. 

Event selection. For further processing, a total of A'tot = 8 895 205 events were selected where at 
least five stations had triggered. Events which fulfilled this condition, but had less than 16 DOMs 
read out (before event cleaning), were reweighted in the analysis with the prescale factor of 5 (see 
above). 

The effective livetime was calculated by fitting the distribution of time differences between 
events. At, with an exponential function, 

N{At) = Nq exp(-AVT). (4) 

This was done individually for each data taking run. The selected runs have a total effective 
hvetime of T X)runs »(^* ' '^^) = (3274.0 ± 1.9) h, which corresponds to 89.4 % of the selected 153 
days period. The uncertainty on the livetime was included in the statistical error. 
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4. Air shower reconstruction 



The energy of the primary particle cannot be measured directly, but has to be determined from 
the air shower parameters. Properties of an air shower that are reconstructed by IceTop are the 
shower core position, its direction, and the shower size. The latter is a measure of primary energy 
and is defined as the signal S^c! measured at a certain distance Ri-^i from the shower axis. These 
properties are reconstructed by fitting the measured charges with a lateral distribution function 
and the signal times with a function describing the geometric shape of the shower front. 

4.I. The reference radius R^^f 

The average logarithmic distance to the shower axis, (logi?), of signals participating in the fit 
for the given array configuration and energy range under investigation is about 125 m. While this 
number does depend on the primary energy and mass, it is limited by the relatively small size and 
the particular geometry of the 26-station array. A constant Rj-d = 125 m was chosen in order to 
minimize the correlation between the parameters S^-cf and /3 in the fit. The shower size parameter 
is thus referred to as 5'i25. 

4-. 2. Time and charge distribution of air shower signals 

Lateral charge distribution. IceTop tanks are not only sensitive to the number of charged particles, 
but also detect photons. Furthermore, the signal generated by a particle when it traverses the tank 
also depends on incident particle type, energy and direction. Therefore, the charge expectation 
value in an IceTop tank at distance R from the shower axis was described by an empirical lateral 
distribution function found in Monte Carlo simulations [9_: 

S{R) = 5,ef • j . (5) 

This is a second order polynomial in logi? for the logarithm of the signal, log S{R): 

log S{R) = log 5rcf - /? log (^^^ - K log2 . (6) 

This function behaves unphysically at small distances to the shower axis {R < Im). However, as 
described in the next subsection, all signals within 11 m of the core, are excluded from the fit. The 
free parameters of the function, in addition to the shower size, S'lof , are /3 and n, corresponding to 
the slope and curvature in the logarithmic representation at i? = -Rrcf- The parameter k, is fixed 
at the average value of 0.303 found in simulation studies and it was verified that this constraint 
does not have a significant impact on the result. Therefore, a fit of function ([6| depends only 
on two explicit parameters (S'rcf, P) and, since R depends on shower core position [x^, tJc) and 
direction {9, 0), implicitly on four more. 

In the following we will only refer to the reference radius of 125 m motivated in the previous 
subsection. Figure |5] shows an example of the lateral distribution function fit of a shower with 25 
triggered stations. 
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Figure 5: Left: Example of an IceTop lateral fit. The shower triggered 25 stations and the reconstructed shower size 
is S'i25 = (65.1 lb 2.8) VEM. Right: Time residuals with respect to a plane perpendicular to the shower direction 
given by Eq. |8]l. "Upstream" and "downstream" refer to tanks being hit before and after the shower core reaches 
the ground. 

Time distribution. The arrival times of the signals map out the shower front. The expected signal 
time of a tank at the position x was thus parametrized as 



t{x) ^to + lixc- x) ■n + At{R). 



(7) 



Here, to is the time the shower core reaches the ground. Xc is the position of the shower core on 
the ground and n is the unit vector in the direction of movement of the shower. The ground was 
defined as the \/S'-weighted average of participating tank altitudes, which varied by about 3 m. 
The term At{R) describes the shape of the shower front as a function of distance R to shower 
axis and is the time residual with respect to a plane perpendicular to the shower axis which 
contains Xc- Experimentally, the shower front can be described by the sum of a parabola and a 
Gaussian function, both symmetric around the shower axis: 

2ct2 



At{R) = a + ( exp ( - ^ ^ 



- 1 



(8) 



with the constants 



4.823 10"^ ns/m^ 



b = -19.41ns, 



a ~ 83.5 m. 



Function ^ is fitted to the measured signal times with five free parameters: two for the core 
position, two for the shower direction and one for the reference time Iq. Hence, the complete air 
shower reconstruction has the following parameters: position of the shower core {xc,yc), shower 
direction 9 and shower size 5*125, slope parameter /3, and time at ground t^. 



4-. 3. Likelihood fit 

Likelihood function. The functions (|6|, (|7| and Q describing the expectations for the charge and 
time of air shower signals were fitted to the measured data using the maximum likelihood method. 
In addition to terms for the signal charges and times, the likelihood function also takes into account 
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stations that did not trigger so that the full likelihood function consisted of three factors. As usual 
we use the logarithm of the likelihood function: 

L=Lq + Ln + Lf (9) 

The first term, 

describes the probability of measuring the charges Si if the fit expectation value at the position of 
the tank is Sf^ as given by the lateral distribution function (|5|. The sum runs over all tanks that 
have triggered. The signal fluctuations are described by a normal distribution of log Si around 
log5f' , with standard deviations aq depending on the signal charge. The charge dependence of aq 
has been determined experimentally from the local shower fluctuations between the two tanks of 
a station and are reasonably well reproduced by simulation [10) . It can roughly be described by 
a linear improvement of log (cTg (log S*)) until a saturation level is reached at S ^ 120 VEM. The 
second sum in Lq accounts for the proper normalization of the signal likelihood and is required 
because the standard deviations aq depend on the fitted signals. 
The next term of the log- likelihood function ([9|, 

Lo = Y.\n[l-{Pf)'), (11) 

j 

accounts for all stations j that did not trigger. The probability that one tank in station j delivers 
a signal at a given charge expectation value is 

' ^T^^oqisf)' J ^"^l 2apfy^ , 

The lower integration limit is defined through the charge threshold of Sj^^' = 0.3 VEM for the tank 
signal, as determined by the retriggering procedure described in Section |3] The charge expectation 
value, S^*", was evaluated for the center of a line joining the centres of the two tanks. Since 
the two tanks of one station are operated in coincidence, there are no single untriggered tanks. 



1 r (logs',- -log s-f) , 



Equation (111 is an approximation because it assumes that in the two tanks is independent. Of 
course, there is a natural correlation in the signal expectation values of two nearby tanks because 
they have a similar value of the lateral distribution function. However, the fiuctuations about this 
expectation value are assumed to be uncorrelated. 

The third term of function ([9|, £t, describes the probability for the measured set of signal 
times, 

= - E T^2(j^ - ^ln(cT,(i?.)/ns), (13) 

where the index i runs over all tanks, ti is the measured signal time of tank i and tf* — t{xi) is 
the fitted expectation value according to function (|7]). The arrival time fiuctuations at{Ri) depend 
on the distance Ri of tank i to the shower axis, and are the RMS of the arrival time distribution 
found in experimental data |10) . 
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Fit procedure. The likelihood fit was seeded with first-guess calculations for the core and the 
direction of the shower. As a first estimate of the core position the centre-of-gravity of tank 
positions Xi weighted with the square root of the charges was calculated: 



The square root of S used as a weight was chosen based on a study of the achievable fit accuracy. 
The starting values for shower direction and arrival time were obtained by fitting a plane to the 
signal times. 

The likelihood minimisation is then done in several iterations to improve the stability of the 
fit. At first the shower direction is fixed and only the lateral fit of the charges is iterated with the 
free parameters 5*125, Pi ^^nd core position. After each iteration, those tanks that are closer than 
11m to the shower axis are removed from the fit. Iteration is stopped when no more tanks are 
removed from the fit. The reason for this step was that very large signals tended to unnaturally 
attract shower cores, which had a negative effect on the shower core resolution in the vicinity of 
stations. Additionally, this mitigated the effect of saturated pulses. Then, a final iteration is done 
in which description of the shower curvature is included and the shower direction is varied. 

5. Simulation of air showers and the IceTop detector 

The relation between the measured signals and the energy of the primary particle, as well as 
detection efficiency and energy resolution were obtained from CORSIKA [11 air shower simulations 
and simulations of the IceTop detector. 

5.1. Air shower simulation 

We simulate the development of air showers in the atmosphere using the simulation code COR- 
SIKA [llj. Inside CORSIKA, the hadronic component of the air showers was simulated using the 
models SIBYLL2.1 [HI [T3] and FLUKA 2008.3 |H [TS] for the high and low energy interactions, 
respectively. The electromagnetic component was simulated using the EGS4 code [16 and no 
'thinning' (reduction of the number of traced particles) was applied. To study systematic effects 
of the hadronic interaction model, small samples of showers were simulated using the QGSJET- 
II |T7l [T5] and EPOS 1.99 [T^ high energy interaction model. Two different parameterizations of 
the South Pole atmosphere from two days in 1997 based on the MSIS-90-E model [SD] were used: 
July 1st and October 1st (CORSIKA atmospheres 12 and 13). The July atmosphere has a total 
overburden of 692.9g/cm^, while the October atmosphere has an overburden of 704.4g/cm^. The 
July atmosphere was used in the data analysis, because its total overburden is close to the average 
measured overburden of 695.5 g/cm^ and its profile corresponds to that of a South Pole winter 
atmosphere. The October atmosphere model was used to study systematic uncertainties due to 
the atmospheric profile used in the simulation. 




(14) 
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5.2. Detector simulation 

The output of the CORSIKA program, i. e. the shower particle types, positions and momenta at 
the observation level of 2835 m, were injected into the Ice Top detector simulation. The simulation 
determines the amount of light produced by the shower particles in the tanks followed by the 
simulation of the PMT, the DOM electronics and the trigger chain. 

The Cherenkov emission inside the tanks is simulated using Geant4 |21l 122]. All structures of 
the tank, the surrounding snow, including individual snow heights on top of each tank, as well 
as the air above the snow are modeled realistically j23j. The snow heights used in the simulation 
corresponded to those measured in January 2007 (see Fig. [4|. In order to save computing time, 
Cherenkov photons are not tracked; only the number of photons emitted in the wavelength interval 
300 nm to 650 nm is recorded. Using Geant4 simulations, that include Cherenkov photon tracking 
until photons reach the PMT, it was shown that the number of detected photons scales linearly 
with the number of emitted photons, independent of incident particle type and energy. The 
propagation of Cherenkov photons is modeled by distributing the arrival times according to an 
exponential distribution, which is tuned such that simulated waveform decay times match those 
observed in experimental data (26.5 ns for zirconium lined tanks and 42.0 ns for tanks with Tyvek 
bag). 

The number of photoelectrons corresponding to 1 VEM was taken from the VEM calibration of 
the real tanks and used as an input for the simulation. The simulated tanks were then calibrated 
by generating muon spectra as in experimental data using air shower simulations with primary 
energies between 3GeV and 30TeV and zenith angles up to 65°. Thus, the ratio between the 
number of emitted Cherenkov photons and observed photoelectrons was determined by the VEM 
calibration of simulated tanks. 

In the next step the generated photoelectrons are injected into a detailed simulation of the 
PMT followed by the analog and digital electronics of the DOM. To simulate the photomultipliers, 
Gaussian single photoelectron waveforms with a random charge according to the average single 
photoelectron spectrum are superimposed [5J. Afterwards, a saturation function is applied to 
the resulting waveforms. In the DOM simulation, the pulse shaping due to the analog front end 
electronics is applied to the output of the PMT simulation. This includes the individual shaping of 
the signal paths to the ATWD and the discriminators, as well as the simulation of the droop effect 
induced by the toroid that couples the high voltage circuits of the PMT to the readout electronics. 
Then, the discriminators are simulated and the local coincidence conditions are evaluated. Finally, 
the waveform digitization and the array trigger are simulated. 

Simulated data are of the same format as the experimental data and were reconstructed in the 
same way, as described in the previous section. 
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5.3. Simulation datasets 

In this analysis we describe the cosmic ray composition just with the two extreme elements 
hydrogen and iron. The justification comes from the fact that the final result is not sensitive to 
details of the composition but only to the mean logarithmic mass. 

In total 2 • lO'^ showers of proton and iron primaries in the energy range between 100 TeV 
and 100 PeV were generated in 30 logarithmic energy bins according to an spectrum. For 
the analysis, the events are reweighted to an E^^ flux, which is closer to the results of previous 



experiments and thus reduces systematic biases (see also Section 8.9 1. In addition to pure proton 
and iron simulations we also combined the datasets using a parametrization of Glasstetter's two- 
component model |[24j . We transformed the proton flux to the form 



dln£; " " I IPeV j I ^ l^knc 



(15) 



as suggested in PS^, with /q = 3.89 • 10"^ ni-^g-isr-i, 71 = -2.67, 72 = -3.39, Skncc = 4.1 PeV, 
and e — 2.1. The iron flux was used as specifled in Ref. [24] : 

ji|£g = l«.10-=„-V.Sr-.(^)"™. (16) 

The total flux was then normalized to the same E~^ spectrum as in case of the single component 
Monte Carlo. 

Since shower generation is CPU intensive the same showers were sampled several times inside 
a circle with a radius of 1200 m around the center of the 26 station IceTop array. The number of 
samples was chosen for different energy bins such that every shower would remain on average only 
once in the flnal sample after applying the cuts described in the next section. This ensures a good 
balance between an effective use of the generated showers and the artiflcial fluctuations introduced 
by oversampling. 



6. Event selection and reconstruction performance 

Quality cuts. Based on the reconstruction results the following quality criteria were required for 
each event entering the final event sample, for both simulated and experimental data: 

• Containment cut: The reconstructed core and the first-guess core position had to be at least 
50 m inside the boundary of the array. The array boundary is defined by the polygon with 
vertices at the centers of stations at the periphery of the array and edges connecting these 
stations. This cut defines a fiducial area of Acut = 0.122 km^. Furthermore, it was required 
that the station containing the largest signal is not on the border of the array. 

• Only events with zenith angles 6 < 46° were considered. 

• The reconstruction uncertainty on the core position had to fulfill (Tcorc = \ r^x ^ ^ 20 m. 
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Figure 6: Relative abundance of proton and iron in our parametrization of Glasstetter's two-component model as a 
function of primary energy. Above about 10 PeV the spectrum is dominated by iron. 



Table 1: Passing rates of the quality cuts described in the text for events with 5i25 > 1 VEM. Statistical errors on 
experimental data are negligible. 



Experimental data Monte Carlo 

Cut 

Passing rate Cumulative Passing rate Cumulative 



A^station > 5 and S'i25 > 1 VEM 


100% 




100% 




Largest signal contained 


42.5% 


42.5% 


(39.4 ±0.5)% 


(39.4 ±0.5)% 


First guess core contained 


95.8% 


40.7% 


(95.4 ±0.4)% 


(37.6 ±0.5)% 


Core contained 


78.9% 


32.1% 


(81.1 ±0.5)% 


(30.5 ±0.6)% 


Zenith 9 < 46° 


96.3% 


30.9% 


(96.4 ±0.5)% 


(29.4 ±0.6)% 


(Tcoro < 20 m 


99.7% 


30.8% 


100% 


(29.4 ±0.6)% 


2.0 < /3 < 4.5 


98.1% 


30.2% 


(99.7 ±0.1)% 


(29.3 ±0.6)% 



• The slope parameter l3 had to be in the range 2.0 < /3 < 4.5 because most events with /3 
values outside this range were badly reconstructed and because /3 was limited in the fit. The 
removed events had predominantly low primary energies, Eq < 1 PeV. 

In the experimental dataset 3 096 334 events passed the quality cuts. Passing rates for the 
individual cuts are shown in Table [T] for events with 5i25 > 1 VEM. Differences between data and 
Monte Carlo are discussed later in Section ISTtI 

Reconstruction performance. Core position and angular resolution, shown in Fig. [7] are key criteria 
for the performance of air shower reconstruction. The Icr core resolution is defined as the 68% 
quantile of the cumulative distribution of the distances between true and reconstructed shower 
cores; correspondingly the angular resolution is defined as the angle between true and reconstructed 
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Figure 7: Core position and angular resolution for showers with 9 < 30° , obtained from the two-component Monte 
Carlo. At high energies, the distance between true and reconstructed core position of 68% of showers is 7 m or less. 
The angle between true and reconstructed directions of 68% of showers at 1 PcV is smaller than 0.8° and this value 
decreases to 0.4° at 100 PeV. 
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Figure 8: Left: Reconstructed shower size spectra in three zenith angle bins. The energy spectrum was derived 
from these spectra using an unfolding method as described in Section[7] On the right, the same spectra are shown, 
weighted with 5^25. In this representation it is clearly visible that the spectra are not pure power laws, but there 
is a clear structure above log5i25 ~ 1-4. 



shower direction. The numbers shown are for showers with zenith angle 9 < 30°, obtained from 
the two-component Monte Carlo after applying the quality cuts listed in the previous paragraph. 
At the highest energies, a core resolution of 7 m and an angular resolution of 0.4° were achieved. 
In the most inclined zenith angle range considere in this analysis, 40° < < 46°, a core resolution 
of 10 m and an angular resolution of 0.5° was achieved. 



7. Determination of energy spectra 

Using the reconstruction methods and quality cuts described in Sections |4] and [6] the shower 
size spectra shown in Fig. [8] were obtained. In this analysis, the data were split into three zenith 
angle ranges roughly equidistant in sec 61, defined as: 

ni = [0°, 30°] , O2 = [30°, 40°] , ^3 = [40°, 46°] . (17) 
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A steepening of the spectral slope is visible at log(S'i25/VEM) = 0.5 and a possible flattening at 
about log(S'i25/VEM) = 1.4. To determine the energy spectrum from measured data, these 5i25 
spectra were unfolded. Unfolding was performed for each zenith angle range independently. 



7.1. General method 

For the unfolding procedure the response of the detector to a primary particle of mass M , 
energy Ep, zenith angle 0, azimuth cf), and core position {xc, Uc) has to be determined from simula- 
tion. In this analysis we consider only an unfolding of energies. Within each zenith angle range, we 
average over the dependencies on zenith, azimuth, and core position. The response of the detector 
is the probability of measuring a shower size 6*125 given a primary energy Ep and mass M in a 
certain zenith range i7fc. 

In a discrete formulation we define a response matrix R which relates the bin contents N.^ [i = 
1, . . . , m) of a measured S'125 spectrum with the bin contents iVJ (j = 1, . . . , n) of a primary energy 
spectrum for a fixed zenith range 0^: 

Nt = iV;. (18) 

ik) 

The response matrix elements R\- ' are defined as acceptance integrals 



E / di?p / dn /dA^ ^M{Ep) (Sl^, I Ep) 

M 

- r,^ r,^ r,. . , r. . " (19) 



^'^ AEi ^cut 



The model fiux $A/(i?p) of nuclei with mass M weighted by their acceptance function 

Pm (^125 I Ep) - p(5i25, nk I Ep, x„ ye, 0; M) (20) 

is integrated over primary energy bin E^, the angles 9 and ip, and area A± projected on a plane 
perpendicular to the particle direction. It is summed over all mass components M that contribute 

(k) 

to the assumed composition model. R^j is normalized to the fiux integrated over bin j in Ep, solid 
angle Qk, and fiducial area Acut- The function pM is the probability of an event with mass M and 
kinematical variables (Ep, x^, yc-, 0, (j)) to be reconstructed with shower size S'125 ™ ^™ * ^^^^ zenith 
angle 9 in the range ilk, and to pass all cuts listed in Section [b] Thus, i?!*^' for a given primary 
energy bin j, is the ratio between number of events measured in S'125 bin i and zenith bin k, that 
pass all cuts, and the true number of events in that energy bin j and zenith bin k inside the 

(k) 

fiducial area. Since the Ep bins of i?^^ are independent, the total fiux model only affects weighting 

(k) 

of events within one bin, but not neighboring bins. The flux normalization in i?^^, cancels out, 



and the dependence on the spectral index of the flux model is small (see also Section 8.9 1. The 



integrals in Eq. ( 19 1 were determined numerically using the Monte Carlo method. 

With the normalisation to the full flux integral the response matrix has the following normali- 
sation properties (we drop the superscript k for zenith range): 



J2R^J^e,, 5]i?.,=l. (21) 
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Figure 9: Response matrix: shower size 5i25 distribution as a function of primary energy for simulated proton (left) 
and iron (right) showers with zenith angles up to 30°. The crosses give the mean value and spread (RMS) of the 
distribution in each energy bin. 



That means, for a given energy bin j the sum of the probabihties to be detected in any signal bin 
is the efficiency e^; for a given 5i25 bin i the probabihty to belong to any energy Ep is unity. The 
efficiency depends on the energies, the core position {xc,yc) and the angles. 



To obtain the primary energy spectrum from the measured signals the matrix equation (18 1 
has to be inverted: 

NJ = iV^ (22) 

For this unfolding procedure we use an iterative algorithm, which properly accounts for the statis- 
tical fluctuations as will be described Section [731 



7.2. Evaluation of response matrices, efficiencies and resolutions 

Figure [9] shows the response matrix for simulated proton and iron primaries in the interval ili 
of smallest zenith angles. In each bin the colour code represents the probability that an event with 
energy Ep yields a signal 5i25. The binning uses a logarithmic scale. 

For computational purposes and to smooth fluctuations in the simulated response matrix, the 
log 5i25 projections of each log Ep bin j were fitted by a normal distribution function yielding the 
mean value (log5i25)j and standard deviation aiogs,j- The normalisation Ej was calculated as the 
ratio between the sum of Monte Carlo event weights in the final sample and the sum of weights of 
events generated inside the fiducial area defined in Section |6j 

^. = ^^- (23) 

Due to migration of shower cores from outside the fiducial area, this quantity can become larger 
than unity. The energy dependence of the parameters (log 6*125), criogSi and e was then fitted by 



empirical functions (see Appendix A I . These functions were used to smooth statistical fluctuations 
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Figure 10: |(a)| Mean shower size as a function of energy for proton showers of various zenith angles. |(b)| Shower size 
ratio between proton and iron showers. The ratio increases for larger zenith angles because the attenuation of iron 
showers is stronger than for proton showers. 
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Figure 1 1 : Total efficiency for proton showers in different zenith angle ranges as a function of energy. The solid lines 
are fits of function ||A.3|| , and the dashed lines are the results of the corresponding fits for pure iron. 
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Figure 12: Energy resolution as defined in the text for proton showers in different zenith angle ranges. The lines 
are fits according to Eq. ||A.2[| in order to guide the eye. 



in the response matrix and to extrapolate the range of simulations to higher energies in order to 
avoid potential artifacts that might be introduced by cutting the spectra off at 100 PeV. 

The mean values and standard deviations are indicated in Fig. [9] by the points with vertical 
bars. The response functions (log 5i25(£'p)) of proton showers for the three different zenith angular 
intervals flk are shown in Fig. |10(a) Since the shower maximum lies above the detector throughout 
the covered energy range, showers from larger zenith angles are more strongly attenuated by the 
atmosphere and thus have a smaller shower size. In Fig. |10(b)[ these points are compared to the 
mean values for iron. The response matrices for proton and iron are very different: on average, 
iron showers have their first interaction at larger height leading to a larger shower age than for 
protons. Iron showers yield a smaller average signal than proton showers with the same primary 
energy. The difference between proton and iron increases at larger zenith angles. This zenith 
angle dependence has been exploited to test the consistency of our data with models for the mass 
composition, as will be discussed in Section |9] 



Figure 11 shows the efficiencies e obtained in the log Ep bins, which are the normalisations of 
the normal distributions of log 5i25 belonging to this bin, for protons and iron nuclei comparing 



all zenith angle intervals. The lines are fits to Eq. (A. 3). Mostly due to the very conservative 



containment criteria, peak efficiencies were significantly below 100%. The maximum efficiencies in 
the three zenith angle ranges flk correspond to the following effective areas: 

: Aeff = (1.051 ± 0.013) • 10^ 
^2 : Acff = (0.900 ± 0.019) • 10^ 
: Aeff = (0.803 ± 0.012) • 10^ m^. 

Within statistical uncertainties the same values were obtained for iron primaries. 



The energy resolution (see Fig. 12 1 has been determined by transforming the log 5125 distribu- 
tion for a given Ep back onto the log i?p axis. It is worst where the detector becomes fully efficient. 
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which happens between 1 and 3 PeV depending on zenith angle. Towards higher energies the res- 
ohition improves, reaching values between 0.04 and 0.12 in log(£'p) at 100 PeV, corresponding to 
a resolution ue/E between 9% and 23%. The improvement of the energy resolution in the thresh- 
old region toward lower energies is a cutoff effect due to the fact that showers of those energies 
will only trigger the detector if they fluctuate upward. This resolution only covers the statistical 
fluctuations, systematic uncertainties are discussed later in Section [8] 

The response matrices obtained by this method depend on the primary composition assumption, 
as well as the hadronic interaction models and the parametrization of the South Pole atmosphere 
assumed in the simulation. 

7.3. Unfolding 

The response matrices were inverted using an iterative unfolding method based on Bayes' 
theorem described in Ref. |26| . which takes into account the total efficiency e and migration due 
to the fluctuations criog(s)- Simply inverting the response matrix R would lead to unnatural 
fluctuations in the result. 

Starting from a prior distribution P]~{Ep^) in the fc-th iteration, the inverse of the response 
matrix R^^ is constructed by inverting P{Si25\Ep^) — R^j using Bayes' theorem: 

Then, an estimate of the energy spectrum Njf. is obtained from the charge spectrum N.f: 

Nlk = ^J2^'P^(4''\Si2,)- (25) 
In the last step of the iteration, P}^{Ep^) is replaced by 

As initial prior, Pq(Ep^) ^ E^^ was chosen. 

After each iteration, the unfolded spectrum was folded with the response matrix, iV/j, = 
RijN^ 1^, and compared to the measured shower size spectrum. A convergence criterion was 

then deflned using the change in between and the measured shower size spectrum 

between two iterations k and fc + 1 , as in |27| : 

Ax^ik, k + l)= x^im,Nn - x'iNk+i,Nn. (27) 
This quantity decreases monotonically during the iteration process. However, at Ax^{k, fc + 1) = 

(k) 

the unfolding would be equivalent to simply inverting i?^ J and the unfolded spectrum would 
fluctuate unnaturally. To avoid this, the iteration was terminated once Ax^(fc,fc + 1) fell below 
a certain value Ax^^j.^^. The value of this limit was determined beforehand using a simple toy 
simulation in which a known spectrum was folded with the response matrix and then, after adding 



22 



statistical fluctuations, unfolded again. In every iteration step of this unfolding procedure, the 
unfolded spectrum was compared to the known true spectrum. Finally, Axterm — 1-1 chosen 
where the agreement with the true spectrum was best on average. 

The error bars on the unfolded spectrum were determined by varying the shower size spectra 
within their statistical errors and repeating the unfolding. This was repeated n — 3000 times and 
the statistical errors in bin j were determined by comparing each unfolding result A^'^t'^' to the 
average result {N'^): 

K)^ = ;^E«^-(^^).)'- (28) 

k=l 

Similarly, bin-to-bin correlations were obtained: 

cov(z,j) ^lY^i^f'^ - iNl^){Nf'^ - (N^),)- (29) 
fc=i 

It was verifled with a simple toy model that this algorithm correctly reproduces a true in- 
put spectrum, which was folded with the detector response and that the error determination is 
correct [25] , 

7.4- Correction for snow 

Snow accumulates constantly on top of the Ice Top tanks, but a manual measurement of the 
snow height is only possible during the austral summer and therefore is done only once every year. 
The detector simulation took the snow depths measured in January, 2007, into account. Data, on 
the other hand, were taken between June and October, 2007, when more snow had accumulated. In 
order to estimate the effect of this difference, the detector response to proton showers with primary 
energies of IPeV, lOPeV and 30PeV and zenith angles 0°, 30° and 40° was simulated assuming 
once the snow heights measured in January 2007 and once those measured in January 2008. In 
January 2007 the average snow depth on top of Ice Top tanks was 20.5 cm, while in January 2008 
the average height on top of the same tanks was 53.2 cm. Assuming constant increase in snow depth 
and proportionality between log 5i25 and snow depth, shower sizes in August, 2007, were estimated. 
This lead to the following zenith angle dependent energy corrections relative to the simulations 
based on the January 2007 snow height measurement, which were applied to all unfolded energy 
spectra. Within the statistical uncertainties, no energy dependence could be observed: 

: A log(i;/PeV) = 0.0368 ± 0.0009, 
^2 : A log(£;/PeV) = 0.0440 ± 0.0013, (30) 
rJa : A log(£;/PeV) = 0.0513 ± 0.0008. 

8. Systematic uncertainties 

All systematic errors are summarized in Table [2] In the following details about the determina- 
tion of the uncertainties of the energy determination and the flux measurement will be given. 
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8.1. Snow height 

To estimate the systematic error due to the energy correction for snow described in Section [7^ 
snow accumulation was assumed proportional to wind speed. The numbers obtained in this way 
were compared to those assuming constant growth of the snow depth (see above). The result of 
this comparison was used as an estimate of the systematic error on energy determination due to 
snow height. 

8.2. Variations of the atmosphere 

As discussed in Section [2^4] variations of the atmosphere affect the observed shower sizes. The 
influence of two parameters of the atmosphere has been studied in a data driven way: the total 
overburden Xq, and the altitude profile dXy{h)/dh. 

First, the days of data taking were ordered according to the total atmospheric overburden Xq. 
Then the 50 days with the highest and the 50 days with the lowest overburden were selected from 
the total of 153 days. The average overburdens during these periods were X\ow — 679g/cm^ and 
-^high — 700g/cm^, yielding a difference of AX ~ 21g/cm^. From the data taken during these 
days shower size spectra were created for each zenith range i7fc. 

By comparing the shower size spectra obtained in the two periods, the dependence of 6*125 with 
atmospheric overburden was derived. The RMS variation of the total atmospheric overburden 
between June 1 and October 31, 2007, of ax^ = 9.86 g/cm^ was used to estimate the systematic 
error on the energy determination due to atmospheric overburden variations. With the given 
statistical precision, an energy dependence of this variation could not be observed. 

In contrast to total overburden the altitude profile of the atmosphere at South Pole undergoes 
a clear annual cycle. To study the effect of varying the atmospheric profile on air shower mea- 
surements the data taking period was divided into a period of very dense atmosphere (July 25th 
to October 10th) and one when the atmosphere was less dense (remaining days between June 1st 
and July 24th and between October 11th and October 31st). Shower size spectra were extracted 
from the data taken in these two periods and by comparing those spectra, an additional systematic 
error due to the atmospheric profile variation was derived. 

8.3. Atmosphere model in simulation 

The CORSIKA simulations used a model of the South Pole atmosphere. A systematic uncer- 
tainty arises from the choice of model since it does not exactly match the average atmosphere 
during the data taking period. To estimate this error on the energy scale simulations two dif- 
ferent atmosphere parametrizations were compared. CORSIKA atmosphere model 12 (July 1, 
1997), which was used in the unfolding procedure, has a total overburden of 692.9 g/cm^ and 
atmosphere model 13 (October 1, 1997) has a total overburden of 704.4g/cm^. Averaging the 
difference in log 6125 between the two simulations above Ep — I PeV the systematic error due to 
the difference of the simulated overburden and the average overburden in data was determined. 



24 



U 



300 
250 
200 
150 
100 
50 



VEM Stability 
Entries 1560 
Mean 1.931e-17 
RMS 0.03056 




-%.5 -0.4 -0.3 -0.2 -0.1 0.1 0.2 0.3 0.4 0.5 

(VEM - (VEM»/<VEM> 



Figure 13: Relative change in the number of photoelectrons corresponding to 1 VEM between consecutive calibration 
runs for all 15 calibration runs and all high and low-gain DOMs. Calibration runs were carried out every two weeks. 
The RMS value of 0.031 of this distribution was used to estimate the systematic error on energy determination. 



8.4- Calibration 

Systematic uncertainties due to calibration can arise for two reasons: variations of the cal- 
ibration constants between calibration runs, and a discrepancy between the calibration of the 
experiment and the detector simulation. 

The first point was addressed by studying the variation of the VEM calibration between cal- 
ibration runs. Figure [13] shows the relative difference in number of photoelectrons corresponding 
to 1 VEM between calibration runs for all DOMs. From the RMS of this distribution the systematic 
uncertainty on the energy reconstruction due to variations of the VEM calibration was estimated 
to be 3.0%. 

The simulated tanks were calibrated using the same procedure as for the real tanks, as described 



in Section 5.2 The conversion factor between Cherenkov photons and photo electrons resulting 
from this calibration has a statistical uncertainty of 1.5%, which was included as a systematic error 
on the energy. 

8. 5. Droop 

The toroid used to decouple the PMT from the signal capture electronics introduces a significant 



droop effect (see Section 2.3 1, which was not corrected for in the analysis. Not correcting for droop 



is not a source of systematic uncertainty in itself if it is done consistently in data and simulation. 
However, discrepancies in the way the droop effect is simulated in the detector Monte Carlo, 
may lead to undesired systematic effects. In order to quantify these effects, the effect of a droop 
correction algorithm on the recorded charges was compared between data and simulation. From 
this comparison a systematic error on the energy determinatino of 1.5% was derived. 
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Figure 14: Shower size ratio for the two-component assumption between SIBYLL and QGSJET and EPOS based 
Monte Carlo simulations for the two component primary assumption and showers with zenith angles up to 30°. 

8.6. PMT Saturation 

Inaccuracies in the simulated saturation behaviour of the PMT could introduce systematic 
uncertainties on the energy determination mostly at high energies. In simulation, saturation sets 
in at higher charges than in the experiment. In order to estimate the effect of this discrepancy on 
the energy spectrum, an artificial, charge-based saturation function was applied to the simulated 
charges to bring the simulated charge spectrum into agreement with experimental data. Then, the 
simulated showers were reprocessed, and the change in log(5i25/VEM) was used to estimate the 
systematic error on the energy. For primary energies below 10 PeV, the systematic error due to 
the difference in saturation behaviour is less than 0.5%. Above 10 PeV it increases exponentially 
to a value of 2.5% at 100 PeV. 

8. 7. Cut efficiencies 

Differences in the effects of quality cuts described in Section [6] when applied to experimental 
and simulated data lead to a systematic uncertainty on the efficiency and consequently on the 
flux normalization. Passing rates of all cuts for data and Monte Carlo events above threshold are 
Hsted in Table [T] There is a relative difference of 3.0% between data and Monte Carlo in the total 
cumulative passing rate, which is included in the systematic uncertainty on the flux. 

8.8. Interaction model 

Small simulation datasets of proton and iron showers created using the high energy hadronic 
interaction models QGSJET-II and EPOS 1.99 in addition to SIBYLL were used to estimate the 
systematic uncertainty due to the modeling of hadronic interactions. Figure [14] shows the shower 
size ratio between SIBYLL and the alternative simulations as a function of primary energy for the 
two-component primary composition assumption and zenith angles up to 30°. Simulations with 
SIBYLL seem to yield systematically smaller shower sizes, and the same observation was made for 
more inclined showers. 
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The systematic error derived in this way is purely based on a comparison of the three interaction 
models. All of these models have different known strengths and weaknesses in their description 
of the underlying physics. Additionally, they all include extrapolations of cross-sections and mul- 
tiplicity distributions to energy ranges not accessible by current collider experiments which are 
relevant in the first few cosmic ray interactions. Thus, there is an unknown systematic error in 
case the range hadronization models does not cover the true behavior. 

8.9. Response matrix 

Limited Monte Carlo statistics introduce uncertainties into the response matrix. Assuming the 
efficiency is constant above the threshold, the flux error induced by uncertainties of the detector 



response can be estimated by the fit error on cq in Eq. (A.3). The uncertainties on the parameters 



ao and 60 in Equations (A.l) and (A. 2) translate to an uncertainty on the energy in the unfolding 
process. These statistical uncertainties on the response matrix were also included in the systematic 
error of the final result. 

Additionally, the flux model used in the simulation also influences the response matrix. A harder 
spectrum leads to larger average shower sizes in an energy bin than a softer one. Simulations based 
on an E^"^ flux and an E^"^ flux were compared with the standard simulation which assumes a 
power law of E^^ . Above the threshold the resulting difference in shower size appears to be 
independent of primary energy. The differences in shower size between the two extreme spectral 
indices were used as an estimate of the systematic error on energy scale due to the assumed flux 
model. 



8.10. Unfolding procedure 

Two parameters besides the response matrix influence the result of the unfolding: the termi- 
nation criterion Ax^^g^^ and the prior distribution Pq- Varying the termination criterion, lead to a 
variation of the total flux, which was included as a systematic error. 

In addition, varying the spectral index of the initial prior Pq between —2.5 and —3.5, a variation 
of the total flux of about 2% was observed. Below the knee region around 3 to 4 PeV, the spectral 
index seems to depend on the prior (in the most inclined zenith interval even up to 10 PeV. Varying 
the prior lead to a variation of the spectral index below the knee in the most vertical zenith band 
by ±0.01, and in the most inclined zenith range by ±0.025. At higher energies variations appear 
to be purely statistical. 

8.11. Summary of systematic errors 

Systematic uncertainties are summarized in Table [2] The total systematic uncertainty was 
determined by quadratically adding the individual contributions. The error on the determination 
of the primary energy in the most vertical zenith angle range is 5.1% below Ep — 10 PeV, and 5.7% 
above. Main contributions are the calibration stability (3.0%), atmosphere (2.7% in total), and 



27 



Table 2: Summary of systematic uncertainties of the energy and flux determination in the three zenith angle intervals 
Oi, fl2, and Q^. The individual points are explained in the text. 



0° < e < 30° 30° < e < 40° 40° < e < 46° 



Uncertainty 


Energy 


r lux 


Energy 


r lux 


Energy 


r lux 


Snow height 


0.4% 




0.4% 




0.4% 




Overburden variation 


0.26% 




1.9% 




3.0% 




Atmosphere profile variation 


2.5% 




1.8% 




1.1% 
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0.9% 
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Droop 


1.5% 




1.5% 




1.5% 




Calibration stability 


3.0% 




3.0% 




3.0% 




Interaction model 


2.1% 




4.3% 




2.0% 




Flux model 


0.7% 




1.0% 




1.0% 




(logs') and (TiogSi25 


0.7% 




1.2% 




0.8% 




Cut passing rates 




3.0% 




3.0% 




3.0% 


Efficiency 




0.9% 




1.6% 




1.2% 


Unfolding procedure 




1.6% 




3.4% 




5.2% 


Total: Ep < 10 PeV 


5.1% 


3.5% 


6.5% 


4.8% 


5.5% 


6.1% 


Ep > 10 PeV 


5.7% 


3.5% 


7.0% 


4.8% 


6.0% 


6.1% 



the hadronic interaction model (2.1%). The systematic influence of the unknown primary compo- 
sition will be discussed in the next section. Furthermore, a flux uncertainty of 3.5% is caused by 
differences in cut efficiencies between data and Monte Carlo, the efficiency calculation in Monte 
Carlo, and the termination criterion and seed in the unfolding procedure. 



9. Energy spectrum 

Figure [15] shows energy spectra for three zenith angular intervals unfolded under three assump- 
tions on the mass composition: all-proton, all-iron and the two-component model |24| explained in 
Section |5.3| The lower end of the energy range of each spectrum was selected where the efficiency 



according to Eq. (A. 3 1 reached 90% of the maximum value. The threshold was determined individ- 
ually for each zenith interval and primary composition assumption. That way the threshold region 
is excluded and the efficiency can be assumed almost constant. Based on the energy resolution, a 
binning of 10 bins per decade was chosen. 

In Fig. |10(b)| it was shown that the difference in shower size between simulated proton and iron 
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Figure 15: Resulting flux measured with IceTop, weighted with E^-"^ . The reconstruction was done using three 
different composition assumptions as described in the text: (a) pure proton, (b) pure iron, and (c) Glasstetter's two- 
component model. In each case, the data were divided into three different zenith angle bands equidistant in sec(6'). 
Based on the assumption of an isotropic flux, the three individual spectra should agree. The boxes indicate the 
systematic errors. 
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Figure 16: All particle spectra obtained with IceTop from air showers with zenith angles up to 30° under three 
different composition assumptions: pure proton, the two-component model, and a mixture of 30% proton and 70% 



showers increases with zenith angle due to the increasing slant depth in the atmosphere, which has 
a different effect for the different masses: iron showers are attenuated more strongly with increasing 
slant depth than proton showers. Since the cosmic-ray flux is isotropic to a few per thousand the 
flux measured in different zenith angular intervals has to be the same. 

In case of the pure proton assumption (Fig. |15(a)[ ) a good agreement between the three spectra 
is observed. Assuming pure iron (Fig. [T5"(b)| , the individual spectra for the three different zenith 
bands clearly disagree at low energies while they start to converge toward higher energies. Agree- 
ment of the three spectra in case of the two-component model (Fig. 15(c)| is good at low and high 
energies. In the intermediate energy range there is some deviation between the spectrum obtained 
from steepest zenith angle range and the other two spectra. However, they are still consistent 
when considering systematic uncertainties. 

Using a comparison of fluxes in each bin of the spectra from the three zenith angle ranges, 
pure iron could be excluded at a >99% confidence level below 25 PeV. This comparison took into 
account both statistical and systemtic errors. The latter were treated in a conservative way by 
assuming no correlations between them for the different zenith angle intervals. Using the same 
comparison and various mixtures of proton and iron, up to 70% of iron cannot be excluded at any 
energy. 

In Fig. [16] the results obtained in the steepest zenith angle range f2i with three primary 
composition assumptions are compared: pure proton, the two-component model, and 70% iron. 
Only the most vertical zenith angle range was chosen, because the difference in size for showers 
initiated by different primaries is smallest in this zenith interval, as seen in Fig. |10(b)[ and because 
systematic uncertainties are smallest in this range. Because the difference in shower size between 
proton and iron decreases toward higher energies, the spectrum obtained under the 70% iron 
assumption is softer than the proton-based result. While the composition model has a sizable 
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influence on the measured all-particle flux below 10 PeV, the difference between the two extreme 
assumptions of pure proton and 70% iron almost disappears above 30 PeV. 

As a flnal result the cosmic ray spectrum is given separately for the assumptions of the pure- 
proton and the two-component model which both yield consistent fluxes in the different zenith angle 
ranges. The systematic errors, as depicted in Fig. [16] by the bands covering the data points, are 
evaluated for all assumptions separately and without including the uncertainty from the unknown 
composition. The 70%-iron case was used in addition to determine the systematic error range on 
the flux due to primary composition. The range of systematic errors lies between the upper border 
of the 70%-iron band and the lower border of the pure-proton band. At 2.4 PeV, for example, the 
allowed fluxes range from 2.65 x IG^^'^ to 3.34 x 10^^"^ GeV^^ m^^ s^^ sr^^. The contribution to 
the systematic uncertainty due to the primary composition decreases from about 30% at 2 PeV to 
less than 1% above about 60 PeV. 



Figure 17 shows the results for pure proton and the two-component model, without the sys- 
tematic error bands, in comparison to a selection of other experiments \29\ l30l \3T\ [32| l33l [34l \35\ 
[MllSZllMllMllinilllllHllS- Table [s] lists the measured fluxes for these two primary composi- 
tion assumptions. The systematic errors on the flux given in the table have been calculated by 
transforming the systematic error on energy into a flux error based on the local spectral index 7: 
AI/I = -fAE/E. This was added quadratically to the systematic error on the flux. 

The two spectra have been fitted with the following parametrization |25] : 



d\x\E 2(T2-Ti)/'^ 

where /knco is the flux at the knee, -Ekncc is the position of the knee, 71 is the spectral index below 
and 72 above the knee, and e describes the sharpness of the knee. In the flt, statistical errors and 



bin-to-bin correlations according to Equations (28) and (29 1 were used. The results are listed in 
Table H 

In pure-proton case, the data points below the knee are not well fltted with the assumption 
of a single slope. This could either be a real feature of the spectrum or an indication of a wrong 
composition because in the region of the flrst two points the energy threshold causes a mass 
dependent efficiency. In order to obtain nevertheless also for the pure-proton case a flt with a 
unique slope below the knee, the flrst two data points have been excluded from the flt. The 
variation of the parameters when including the flrst or the second point respectively was used as 
a systematic error. When including all data points, all parameters lie within this range, but only 
a bad fit is achieved. 

The systematic uncertainty of the knee energy -Eknco is the systematic error on energy determi- 
nation at that primary energy as given in Section [8] The systematic error of /knee has been obtained 
by quadratically adding the systematic error on the flux determination and the systematic energy 
error transformed into a flux uncertainty based on the local spectral index. In order to determine 
systematic errors on 71 and 72, the fit was repeated using the systematic errors of the data points 
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Table 3: All-particle cosmic ray energy spectra measured by the IceTop air shower array for the pure proton and 
the two-component primary composition assumptions using the hadronic interaction model SIBYLL2.1. 



Energy dN/dE ± stat ± syst (GeV"^ m'^ s"^ sr-^) 



(10® GeV) Proton Two-component 
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(2.098 ±0.023 ±0.4) 


X 


10- 


-14 


7.71 


(0.910 ±0.017± 0.15) 


X 


10- 


-14 


(1.021 ±0.015 


±0.20) 


X 


10- 


-14 


9.70 


(4.41 


±0.10 


±0.8) 


X 


10- 


-15 


(4.92 


±0.10 


±0.9) 


X 


10- 


-15 


12.21 


(2.13 


±0.07 


±0.4) 


X 


10- 


-15 


(2.38 


±0.06 


±0.4) 


X 


10- 


-15 


15.38 


(1.08 


±0.04 


±0.18) 


X 


10- 


-15 


(1.177 ±0.04 


±0.20) 


X 


10- 


-15 


19.36 


(5.01 


±0.26 


±0.9) 


X 


10- 


-16 


(5.58 


±0.23 


±1.0) 


X 


10- 


-16 


24.37 


(2.45 


±0.17 


±0.4) 


X 


10- 


-16 


(2.66 


±0.15 


±0.5) 


X 


10- 


-16 


30.68 


(1.44 


±0.11 


±0.22) 


X 


10- 


-16 


(1.51 


±0.10 


±0.23) 


X 


10- 


-16 


38.62 


(7.0 


±0.7 


±1.2) 


X 


10- 


-17 


(7.5 


±0.7 


±1.3) 


X 


10- 


-17 


48.62 


(3.6 


±0.5 


±0.6) 


X 


10- 


-17 


(3.72 


±0.4 


±0.6) 


X 


10- 


-17 


61.21 


(1.91 


±0.29 


±0.29) 


X 


10- 


-17 


(1.97 


±0.25 


±0.3) 


X 


10- 


-17 


77.06 


(1.04 


±0.18 


±0.19) 


X 


10- 


-17 


(1.05 


±0.17 


±0.19) 


X 


10- 


-17 


97.01 


(4.6 


±1.1 


±1.0) 


X 


10- 


-18 


(4.7 


±1.0 


±1.0) 


X 


10- 


-18 



32 



> 



X 



10" 




log(E/PeV) 



Figure 17: The all-particle cosmic ray energy spectrum obtained from the analysis of IceTop data of events with 
zenith angles up to 30° compared to a selection of other experimental results 1291 1301 1311 1321 1331 1341 1351 1361 1371 1381 
[39l|l0l|lll|42l|43l. 

Table 4: Fit parameters of the cosmic-ray energy spectrum according to function ( |31| for the pure proton and 
two-component model primary composition assumptions. Systematic errors were derived as described in the text 
and exclude the systematic error due to the unknown composition, since these are fits of spectra derived under 
specific composition assumptions. 

(a) Proton 

Parameter Best fit 

4„ec/10-"m-2s-isr-i 

Skncc/PeV 

71 

72 

e 



3.8 ± 1.9(stat) ii.3(syst) 

3.2 ±0.9(stat) to.2(syst) 

-2.5 ± 0.4(stat) j:°;^(syst) 

-3.076 ±0.019(stat) ± 0.15(syst) 
6 ± 4(stat) 
15.6/12 



(b) Two Components 



Parameter 



Best fit 



/k„cc/10-^m-2s-isr-i 

Skncc/PeV 

7i 

72 

e 



2.38 ± 0.23(stat) ± 0.5(syst) 

4.32 ±0.22(stat) ± 0.18(syst) 

-2.759 ±0.015(stat) ± 0.21(syst) 

-3.107 ± 0.016(stat) ± 0.3(syst) 

9 ± 3(stat) 
19.4/13 
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Table 5: Parameters of the flattening of the spectrum at high energy. The errors given are only statistical. 



Parameter Proton Two Components 

-Bbreak/PeV 21 ± 4 23 ± 5 

73 -2.82 ±0.10 -2.87 ±0.09 

xV^df 6.1/10 7.1/11 



as statistical errors. In case of the proton assumption the systematic uncertainty introduced by 
the removal of the first two data points from the fit (see above) has been added quadratically to 
these numbers, which increases the systematic values for this assumption, in particular those of 
the slope 71 below the knee and the knee energy. 

Above about 22 PeV a possible fiattening of the spectrum can be observed independent of pri- 
mary composition assumption. This feature is also visible in the measured shower size spectra 
(see Fig. [Sj. In order to test its statistical significance, the spectra were fitted with function (311 
plus an additional hard break at inbreak with spectral index 73. The goodness of fit improves 
to /Ndf = 6.1/10 for the pure-proton assumption and to X^/^di = 7.1/11 for the two-component 
model assumption. These improvements of the correspond to significances of 2.7 and 3.2 stan- 
dard deviations respectively, which, however, does not include systematic errors. The parameters 
of the flattening are listed in Table |5] 



10. Summary 

We have derived the all-particle cosmic ray energy spectrum in the energy range between 1 PeV 
and 100 PeV from data taken between June and October 2007 with the 26-station configuration of 
the IceTop air shower array at South Pole. 

Using the air shower simulation package CORSIKA with the high-energy hadronic interaction 
model SIBYLL2.1 the relation between shower size 5*125 and primary energy, as well as the detection 
efficiency and energy resolution were determined. Three different assumption on the primary mass 
composition were used as input: pure proton, pure iron and a simple two-component model |24| . 
Based on these results, shower size spectra obtained in three zenith angle ranges were unfolded 
with a Bayesian unfolding algorithm to obtain energy spectra. 

In case of pure proton and the two-component model, it was found that the spectra obtained in 
the different zenith angle ranges were in good agreement. In the pure iron case, on the other hand, 
a strong disagreement between the three spectra was observed at low energies. Since one can safely 
assume that cosmic rays are isotropic in the given energy range, the spectra in all three zenith 
angle ranges should be the same. With this assumption, we concluded that pure iron primaries 
can be excluded below energies of 25 PeV. 

We showed that the attenuation of air showers with increasing zenith angle bears exploitable 
information about the chemical composition of cosmic rays. Nevertheless, the main source of 
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systematic error still remains the primary mass composition. The systematic error due to the 
choice of a hadronic interaction model is relatively small in this analysis because most air shower 
signals are dominated by the electromagnetic component of an air shower, which is relatively well 
understood. For the final result, only the spectra obtained from the most vertical zenith angle 
range, 0° < 6* < 30°, were considered because in this range the dependence on composition and 
systematic errors are smallest. 

In case of the pure-proton assumption the knee in the cosmic-ray energ spectrum was observed 
at 3.2 PeV with a spectral index of —2.5 below and —3.08 above the knee. For the two-component 
model assumption the knee position was determined at 4.3 PeV with spectral indices of —2.76 
below and —3.11 above. Around an energy of 22 PeV an indication of a flattening of the cosmic 
ray spectrum to an index of about —2.85 was observed at the 3(t level. 

Since the completion of IceTop and IceCube in 2011, the array is three times larger than the 
configuration used in this analysis. With this larger array, statistics and containment of high-energy 
showers will be much better, allowing to extend the analysis to higher energies. The main strength 
of IceTop, however, is the possibility to measure air showers at the surface in coincidence with 
high-energy muons penetrating deep enough into the ice to trigger IceCube. The ratio between 
the two measurements is highly sensitive to the mass of the primary particle. 



Appendix A. Parametrization of the response matrix 

In order to mitigate the effects of statistical fluctuations in the unfolding procedure, the response 



matrices described in Section 7.2 were separated into mean logarithmic shower size (logS'125), 
resolution ciogSi ^^'^ efficiency e. There dependences on x — log Ep were then fitted by empirical 
functions: 

(logS'125) (a;) ^ ao + x+ 

I exp{aix) + exp(a2 + asx + a4X^)\ 

In 1 . . ^ , (A.l) 

\^ 1 + exp(a2) J 

60(1 +cxp(&364)) -hexp(-fei)(exp(-62a;) - 1) , , 

1 -hexp(-&3(x - bi)) 

Co 

J X < C4 

1 + exp(-ci(a; - C2) -|- C3(x - €4)2) 
e{x) = { . (A. 3) 

Co . 

X > C4 



and 



1 + exp(-ci(a; - C2)) 
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